Comparative Genomics Southern Flounder

knitr::opts_chunk$set(warning = FALSE, message = FALSE)

source("lib/libraries.R")
source("lib/ggplot.R")

Family A Female Map

Load data

Required inputfiles are genfile: *.raw generated using convert.pl to convert filtered *.tsv-output from haplotyper and mapfile: *.map text file with two columns (first column is putative linkage group, second column is markername). Need to create mapfile based on markers in genfile before analysis.

# load data
FamA_F <- read.cross(format = "mm", 
                     dir = "data/LINKMAP/", 
                     file = "FamA-F.raw", 
                     mapfile = "FamA-F.map", 
                     estimate.map = FALSE)

# number of individuals/loci in data set
Sum_FamA_F.raw <- summary(FamA_F)

Raw data set contains 162 individuals and 1383 loci.

Pre-mapping quality control

Filter 1: Missing data

Missing genotypes (loci/individuals)

Omit loci with > 5% missing data and individuals with > 25% missing data.

After filtering data set contains 162 individuals and 1364 loci.

Filter 2: Omit duplicate individuals

Compare genotypes of pairs of individuals to identify pairs with unusually similar genotypes

No matching individuals were identified.

# # identify pairs with over 90% matching genotypes
# match.i <- which(cg > 0.9, arr.ind=TRUE)
# match.i <- match.i[match.i[,1] < match.i[,2],] # table with row (genotype 1) and column (genotype 2) that were found to match in cg
# match.i 
# if no matching individuals identified
FamA_F.f2 <- FamA_F.f1
# number of individuals/loci in data set
Sum_FamA_F.f2 <- summary(FamA_F.f2)
Nind_F.f2 <- as.numeric(Sum_FamA_F.f2$n.ind)
Nloci_F.f2 <- sum(as.numeric(Sum_FamA_F.f2$n.mar))

Data set contains 152 individuals and 1468 loci.

Filter 3: Omit duplicate markers

Identify matching markers. Drop duplicate markers but create table of corresponding marker to be mapped so the can be added post-mapping.

# identify matching markers (matching genotypes/including missing data)
match.loci <- findDupMarkers(FamA_F.f2, exact.only = TRUE)
# keep table to add markers back in after mapping
match.loci_F <- ldply(match.loci, data.frame) %>%
  rename(LOCUS = `.id`, MATCHLOC = `X..i..`)
write.table(match.loci_F, file = "data/LINKMAP/MatchLociFamAF.txt", 
            quote = FALSE, row.names = FALSE, sep = "\t")
# create list of matching loci
match.l <- as.character(match.loci_F$MATCHLOC) 
# drop duplicate markers
FamA_F.f3 <- drop.markers(FamA_F.f2, match.l)
# number of individuals/loci in data set
Sum_FamA_F.f3 <- summary(FamA_F.f3)
Nind_F.f3 <- as.numeric(Sum_FamA_F.f3$n.ind)
Nloci_F.f3 <- sum(as.numeric(Sum_FamA_F.f3$n.mar))

After filtering 1246 loci are retained.

Filter 4: Drop loci with distroted segregation patterns

Loci should have specific segregation patterns based on marker type.

# calculate genotype frequencies and p-value (after Bonferroni correction)
geno.freq.l <- geno.table(FamA_F.f3)
# create list of markers with p-value <0.05
seg.dist.l <- geno.freq.l[geno.freq.l$P.value < 0.05/totmar(FamA_F.f3),] 
drop.l.seg.dist <- rownames(geno.freq.l[geno.freq.l$P.value < 1e-10,])
# drop markers with significant (p<0.05) deviation from expected segregation patterns
FamA_F.f4 <- drop.markers(FamA_F.f3, drop.l.seg.dist)
# number of individuals/loci in data set
Sum_FamA_F.f4 <- summary(FamA_F.f4)
Nind_F.f4 <- as.numeric(Sum_FamA_F.f4$n.ind)
Nloci_F.f4 <- sum(as.numeric(Sum_FamA_F.f4$n.mar))

After filtering 1241 loci are retained.

Filter 5: Genotype frequencies by individual

Determine genotype frequencies by individual and remove individuals as necessary.

After filtering 152 individuals are retained.

Determine correct linkage phase

Estimate recombination fraction (rf) & LOD score for each marker pair. Potentially switched markers are pairs of marker with strong association (LOD score) but rf >> 1/4.

Form initial linkage groups and re-organize loci into initial LGs.

# form initial linkage groups
init.lg <- formLinkageGroups(FamA_F.LP, max.rf = 0.35, min.lod = 6)
table(init.lg[,2])

 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 
36 30 28 27 26 26 25 24 23 23 23 23 23 22 21 20 20 20 20 20 20 20 20 19 19 19 19 19 
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 
18 18 18 18 18 17 17 17 16 16 15 15 15 13 13 11 11 11 11  7  1 

# re-organize markers into inferred initial LGs
FamA_F.LP <- formLinkageGroups(FamA_F.LP, max.rf = 0.35, min.lod = 6, reorgMarkers=TRUE)

Plot rf vs LOD to identify chromosomes to switch alleles for.

Switch alleles, re-estimate rf, form new linkage groups and re-organize markers.

# use plot to identify chromosomes to switch alleles for
for (chr in c(4:5, 10, 18, 21:23, 28:29,
              31:32, 34:37, 39:41, 43, 45:48)){
  toswitch <- markernames(FamA_F.LP, chr=(chr))
  FamA_F.LP <- switchAlleles(FamA_F.LP, toswitch)
}

# re-estimate rf
FamA_F.LP <- est.rf(FamA_F.LP)

# form new linkage groups & re-organize markers
init.lg <- formLinkageGroups(FamA_F.LP, max.rf=0.35, min.lod=6)
table(init.lg[,2])
FamA_F.LP <- formLinkageGroups(FamA_F.LP, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)

Plot rf vs LOD:

Switch alleles, re-estimate rf, form new linkage groups and re-organize markers.

for (chr in c(23)){
  toswitch <- markernames(FamA_F.LP, chr=(chr))
  FamA_F.LP <- switchAlleles(FamA_F.LP, toswitch)
}
# re-estimate rf
FamA_F.LP <- est.rf(FamA_F.LP)

Alleles potentially switched at markers 
  dDocent_Contig_11546

# form new linkage groups & re-organize markers
init.lg <- formLinkageGroups(FamA_F.LP, max.rf=0.35, min.lod=6)
table(init.lg[,2])

 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
62 55 48 48 46 44 43 42 41 41 40 37 37 36 36 35 35 35 34 34 32 27 22 20  1 

FamA_F.LP <- formLinkageGroups(FamA_F.LP, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)

Plot rf vs LOD:

Drop marker in LG25 and 26:

Establish initial marker order

Establish initial marker order within LGs using greedy algorithm. Adds one marker at a time with previously added marker locations fixed. Position determined by minimizing number of obligate cross overs.

# Create final cross/map object
FamA_F.final <- FamA_F.LP

# establish initial marker orders
FamA_F.final <- orderMarkers(FamA_F.final, chr = 1:24)

Initial quality assessment

Examine pairwise linkage & Identify gaps between markers

Plot rf vs LOD after ordering markers:

Large gaps between markers indicate that adjacent markers are only weakly linked. Compare average and largest spacing per LG.

summaryMap(FamA_F.final)

        n.mar length ave.spacing max.spacing
1          62   63.4         1.0         5.9
2          55  107.4         2.0        46.4
3          48   71.4         1.5         6.4
4          48  123.2         2.6        65.2
5          46   54.2         1.2         6.6
6          44   63.7         1.5        12.5
7          43   56.0         1.3         6.6
8          42   60.0         1.5         6.6
9          41   49.6         1.2         8.8
10         41   86.6         2.2        18.4
11         40   63.9         1.6        14.2
12         37   65.7         1.8        10.2
13         37   49.6         1.4         6.6
14         36   61.0         1.7        10.3
15         36   57.7         1.6         8.7
16         35   52.0         1.5         5.9
17         35   53.4         1.6        11.0
18         35   91.4         2.7        23.1
19         34   84.8         2.6        31.6
20         34   55.1         1.7         7.3
21         32   46.0         1.5         5.2
22         27   58.4         2.2        16.7
23         22   57.6         2.7        16.7
24         20   35.3         1.9        22.2
overall   930 1567.4         1.7        65.2

Drop markers one by one to determine change in chromosome length & change in log likelihood. Identify markers that result in smaller estimated chromosome length and increased LOD. Dropping terminal markers will usually result in decreased length.

# drop each marker and compare effect
dropone <- droponemarker(FamA_F.final, maxit = 4000, tol = 1e-6, 
                         error.prob = 0.001, verbose = TRUE) %>%
  add_rownames("Marker") %>%
  gather("Effect", "Change", 4:5)

# plot results
plot_dropone <- function(dropone){
  p <- ggplot(dropone, aes(x = pos, y = Change)) +
    geom_line() +
    geom_point(shape = 21, size = 2,
               color = "black", fill = "white") +
    facet_grid(Effect ~ ., scales = "free") +
    labs(x = paste("LG", dropone$chr, "Position", sep =" "), y = "Change if marker dropped") +
    theme_facet
  print(p)
}

dropone %>%
  group_by(chr) %>%
  do(plot = plot_dropone(.)) %>%
  ungroup()

Compare effect in terms of change in length and likelihood and identify top 10 loci with strongest effect per LG.

# data frame with top10 markers with strongest effect
dropone_sum <- dropone %>%
  spread("Effect", "Change") %>%
  group_by(chr) %>%
  top_n(10, LOD) %>%
  ungroup() %>%
  arrange(chr, LOD)

Drop markers as necessary and re-estimate map.

# drop markers (LOD more positive, larger diff in likelihood)
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_50993") # LG24
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_49388") # LG19
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_38538") # LG14
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_37810") # LG14
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_27419") # LG12
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_328929") # LG10
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_42202") # LG9
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_103540") # LG4
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_85985") # LG4
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_3449") # LG4
FamA_F.final <- drop.markers(FamA_F.final, "dDocent_Contig_5061") # LG2

# re-estimate map after eliminating markers
FamA_F.dropl <- est.map(FamA_F.final, error.prob=0.001)
FamA_F.final <- replace.map(FamA_F.final, FamA_F.dropl)

# plot maps of all chromosomes
summaryMap(FamA_F.final)
par(ps=12, cex.lab=1)
plotMap(FamA_F.final, horizontal=FALSE, shift=TRUE, show.marker.names=FALSE, alternate.chrid=FALSE)

Identify individuals with high number of XO

Observed number of crossovers per individual:

# remove individuals with excessively high number of crossovers (< 50)
FamA_F.final <- subset(FamA_F.final, ind=(countXO(FamA_F.final) < 50))

Estimate genotyping error rate

Estimate genotyping error rate.

# Estimate genotyping error rate
loglik <- err <- c(0.0001, 0.0005, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.015, 0.02)
for(i in seq(along=err)) {
  cat(i, "of", length(err), "\n")
  tempmap <- est.map(FamA_F.final, error.prob=err[i])
  loglik[i] <- sum(sapply(tempmap, attr, "loglik"))
}

1 of 9 
2 of 9 
3 of 9 
4 of 9 
5 of 9 
6 of 9 
7 of 9 
8 of 9 
9 of 9 

lod <- (loglik - max(loglik))/log(10)
# plot log10 likelihood
par(mfrow=c(1,1), ps=12, cex=1)
plot(err, lod, xlab="Genotyping error rate", xlim=c(0,0.02), ylab=expression(paste(log[10], "likelihood")))

Genotyping error rate with highest likelihood is 0.0005.

Identify genotyping errors and code as missing

Tight double-crossovers (single marker out of phase with adjacent markers) indicate genotyping errors. Identify double cross-overs by calculating genotyping error LOD scores (Lincoln & Lander 1992). Calculate errorload (genotype being in error vs not being in error). Plot grid of LOD scores indicating genotypes likely to be in error: darker pixels likely to be error: white (LOD <2), grey (2 < LOD < 3), pink (3 < LOD < 4.5), purple (LOD >4.5)

Eliminate genotypes (coded as missing) for LOD > 4 and re-estimate map.

# determine cut-off (LOD > 6)
toperr <- top.errorlod(FamA_F.final, cutoff = 4)

write.table(toperr, "data/LINKMAP/FamA_F.xo",
            col.names = TRUE, quote = FALSE)

# look at genotypes flagged on selected chromosomes
# par(mfrow[1,1], ps=12, cex=1)
# plotGeno(FamA_F.final, chr=1, ind=toperr$id[toperr$chr==1], cutoff=6, include.xo=FALSE)

# eliminate genotypes if necessary (make them missing)
FamA_F.final <- FamA_F.final
for(i in 1:nrow(toperr)) {
  chr <- toperr$chr[i]
  id <- toperr$id[i]
  mar <- toperr$marker[i]
  FamA_F.final$geno[[chr]]$data[FamA_F$pheno$id==id, mar] <- NA
}

# re-estimate map
FamA_F.dropl <- est.map(FamA_F.final, error.prob = 0.001)
FamA_F.final <- replace.map(FamA_F.final, FamA_F.dropl)

247 genotypes were coded as missing.

Re-order markers by LG and compare alternate orders

Compare number obligate XO, LOD & Chromosome length to determine best order

Re-order markers after loci and individuals have been removed and use sliding window to compare alternate orders in terms of number of obligate crossovers and likelihood.

# compare number of obligate crossovers for alternate orders
FamA_F.XO <- lapply(c(1:24), function(chr) {
  rip <- ripple(FamA_F.final, chr = chr, window = 7, method="countxo", error.prob = 0.001, maxit = 2000, tol = 1e-6)
  rip <- rip[1:11,]
})

# access each element in list using [[]]
# Ch1_XO <- as.data.frame(FamA_F.XO[[1]]) %>%
#  rename(Chr1 = obligXO) %>%
#  select(Chr1)

FamA_F.lik <- lapply(c(1:24), function(chr) {
  rip.lik <- ripple(FamA_F.final, chr = chr, window = 5, method="likelihood", error.prob=0.002, maxit=2000, tol=1e-4)
  rip.lik <- rip.lik[1:11,]
})

Switch marker order as necessary

Identify order with lowest number of XO and highest LOD while minimizing chromosome length.

change.marker.order <- function(LG, neworder) {
  Orderlikelih <- as.data.frame(FamA_F.lik[[LG]]) %>%
    select(-LOD) %>%
    select(-chrlen)
  switch_marker <- as.numeric(Orderlikelih[3,])
  FamA_F.final <- switch.order(FamA_F.final, chr = LG, switch_marker, error.prob=0.001)
}

## LG1
# Identify order with lowest no. of XO & highest LOD while minimizing chromosome length
FamA_F.XO[[1]]
FamA_F.lik[[1]]
# change.marker.order(1, new)

## LG2 
FamA_F.XO[[2]]
FamA_F.lik[[2]]
change.marker.order(2, 1)

## LG3
FamA_F.XO[[3]]
FamA_F.lik[[3]]
# change.marker.order(3, 1)

## LG4
FamA_F.XO[[4]]
FamA_F.lik[[4]]
# change.marker.order(4, neworder)

## LG5
FamA_F.XO[[5]]
FamA_F.lik[[5]]
# change.marker.order(5, 2)

## LG6
FamA_F.XO[[6]]
FamA_F.lik[[6]]
# change.marker.order(6, neworder)

## LG7 
FamA_F.XO[[7]]
FamA_F.lik[[7]]
# change.marker.order(7, neworder)

## LG8
FamA_F.XO[[8]]
FamA_F.lik[[8]]
# change.marker.order(8, neworder)

## LG9
FamA_F.XO[[9]]
FamA_F.lik[[9]]
# change.marker.order(9, neworder)

## LG10
FamA_F.XO[[10]]
FamA_F.lik[[10]]
change.marker.order(10, 1)

## LG11
FamA_F.XO[[11]]
FamA_F.lik[[11]]
# change.marker.order(11, new)

## LG12
FamA_F.XO[[12]]
FamA_F.lik[[12]]
# change.marker.order(12, 1)

## LG13
FamA_F.XO[[13]]
FamA_F.lik[[13]]
# change.marker.order(13, new)

## LG14
FamA_F.XO[[14]]
FamA_F.lik[[14]]
# change.marker.order(14, new)

## LG15
FamA_F.XO[[15]]
FamA_F.lik[[15]]
change.marker.order(15, 1)

## LG16
FamA_F.XO[[16]]
FamA_F.lik[[16]]
# change.marker.order(16, neworder)

## LG17
FamA_F.XO[[17]]
FamA_F.lik[[17]]
# change.marker.order(17, neworder)

## LG18
FamA_F.XO[[18]]
FamA_F.lik[[18]]
# change.marker.order(18, 1)

## LG19
FamA_F.XO[[19]]
FamA_F.lik[[19]]
# change.marker.order(19, new)

## LG20
FamA_F.XO[[20]]
FamA_F.lik[[20]]
# change.marker.order(20, new)

## LG21
FamA_F.XO[[21]]
FamA_F.lik[[21]]
change.marker.order(21, 4)

## LG22
FamA_F.XO[[22]]
FamA_F.lik[[22]]
# change.marker.order(22, new)

## LG23
FamA_F.XO[[23]]
FamA_F.lik[[23]]
# change.marker.order(23, new)

## LG24
FamA_F.XO[[24]]
FamA_F.lik[[24]]
# change.marker.order(24, new)

# re-estimate map after re-ordering all LGs (chromosomes)
FamA_F.reorder <- est.map(FamA_F.final, error.prob = 0.001)

# re-place map portion of cross
FamA_F.reorder <- replace.map(FamA_F.final, FamA_F.reorder)

Post-odering quality assessment II

Examine pairwise linkage & Identify gaps between markers

Plot rf vs LOD

Identify large gaps between markers which is indicative of adjacent markers being only weakly linked.

summaryMap(FamA_F.final)

        n.mar length ave.spacing max.spacing
1          62   61.0         1.0         5.9
2          54  106.5         2.0        46.6
3          48   62.8         1.3         5.9
4          45   56.5         1.3         7.3
5          46   54.1         1.2         6.6
6          44   61.2         1.4        11.5
7          43   53.9         1.3         6.6
8          42   59.3         1.4         6.6
9          40   47.3         1.2         8.8
10         40   59.0         1.5         7.3
11         40   62.8         1.6        14.2
12         36   54.6         1.6        10.3
13         37   48.8         1.4         6.6
14         34   40.0         1.2         7.3
15         36   56.2         1.6         8.8
16         35   51.5         1.5         5.9
17         35   52.7         1.6        11.0
18         35   90.3         2.7        23.2
19         33   51.6         1.6         6.6
20         34   54.9         1.7         7.3
21         32   45.9         1.5         5.2
22         27   56.9         2.2        16.4
23         22   57.5         2.7        16.7
24         19   11.5         0.6         4.5
overall   919 1356.8         1.5        46.6

Plot maps of all linkage groups:

Drop markers one by one determine change in chromosome length & change in log likelihood. Identify markers that result in smaller estimated chromosome length and increased LOD.

# drop each marker and compare effect
dropone <- droponemarker(FamA_F.final, maxit = 4000, tol = 1e-6, 
                         error.prob = 0.001, verbose = TRUE) %>%
  add_rownames("Marker") %>%
  gather("Effect", "Change", 4:5)

# plot results
plot_dropone <- function(dropone){
  p <- ggplot(dropone, aes(x = pos, y = Change)) +
    geom_line() +
    geom_point(shape = 21, size = 2,
               color = "black", fill = "white") +
    facet_grid(Effect ~ ., scales = "free") +
    labs(x = paste("LG", dropone$chr, "Position", sep =" "), y = "Change if marker dropped") +
    theme_facet
  print(p)
}

dropone %>%
  group_by(chr) %>%
  do(plot = plot_dropone(.)) %>%
  ungroup()

Compare effect in terms of change in length and likelihood and identify top 10 loci with strongest effect per LG.

# data frame with top10 markers with strongest effect
dropone_sum <- dropone %>%
  spread("Effect", "Change") %>%
  group_by(chr) %>%
  top_n(10, LOD) %>%
  ungroup() %>%
  arrange(chr, LOD)

Plot rf vs LOD

Identify individuals with high number of XO

Finalize map

# print summary map and map data
SummaryFamA_F <- summaryMap(FamA_F.final)

write.table(SummaryFamA_F, file = "results/SummaryFamA_F.map",
          quote = FALSE, sep = "\t", row.names = FALSE)

FamA_F.map <- pull.map(FamA_F.final,as.table = TRUE)

write.table(FamA_F.map, file = "data/LINKMAP/FamA_F.map",
            quote = FALSE, sep = "\t", row.names = TRUE)

Write files of finalized maps. Write files of finalized map. Map contains 1321 loci grouped into 24 LG.

Family A Male Map

Load data

Required inputfiles are genfile: *.raw generated using convert.pl to convert filtered *.tsv-output from haplotyper and mapfile: *.map text file with two columns (first column is putative linkage group, second column is markername). Need to create mapfile based on markers in genfile before analysis.

# load data
FamA_M <- read.cross("mm", dir="data/LINKMAP/", 
                     "FamA-M.raw", "FamA-M.map", estimate.map = FALSE)

 --Read the following data:
    Type of cross:          bc 
    Number of individuals:  162 
    Number of markers:      1379 
    Number of phenotypes:   1 
 --Cross type: bc 

# number of individuals/loci in data set
Sum_FamA_M.raw <- summary(FamA_M)

Raw data set contains 162 individuals and 1379 loci.

Pre-mapping quality control

Filter 1: Missing data

Missing genotypes (loci/individuals)

Omit loci with > 5% missing data and individuals with > 25% missing data.

After filtering data set contains 162 individuals and 1364 loci.

Filter 2: Omit duplicate individuals

Compare genotypes of pairs of individuals to identify pairs with unusually similar genotypes

No matching individuals were identified.

# # identify pairs with over 90% matching genotypes
# match.i <- which(cg > 0.9, arr.ind=TRUE)
# match.i <- match.i[match.i[,1] < match.i[,2],] # table with row (genotype 1) and column (genotype 2) that were found to match in cg
# match.i 

# if no matching individuals identified
FamA_M.f2 <- FamA_M.f1

# number of individuals/loci in data set
Sum_FamA_M.f2 <- summary(FamA_M.f2)
Nind_M.f2 <- as.numeric(Sum_FamA_M.f2$n.ind)
Nloci_M.f2 <- sum(as.numeric(Sum_FamA_M.f2$n.mar))

Data set contains 152 individuals and 1468 loci.

Filter 3: Omit duplicate markers

Identify matching markers. Drop duplicate markers but create table of corresponding marker to be mapped so the can be added post-mapping.

# identify matching markers (matching genotypes/including missing data)
match.loci <- findDupMarkers(FamA_M.f2, exact.only = TRUE)

# keep table to add markers back in after mapping
match.loci_M <- ldply(match.loci, data.frame) %>%
  rename(LOCUS = `.id`, MATCHLOC = `X..i..`)

write.table(match.loci_M, file = "data/LINKMAP/MatchLociFamAM.txt", 
            quote = FALSE, row.names = FALSE, sep = "\t")

# create list of matching loci
match.l <- as.character(match.loci_M$MATCHLOC) 

# drop duplicate markers
FamA_M.f3 <- drop.markers(FamA_M.f2, match.l)

# number of individuals/loci in data set
Sum_FamA_M.f3 <- summary(FamA_M.f3)
Nind_M.f3 <- as.numeric(Sum_FamA_M.f3$n.ind)
Nloci_M.f3 <- sum(as.numeric(Sum_FamA_M.f3$n.mar))

After filtering 1249 loci are retained.

Filter 4: Drop loci with distroted segregation patterns

Loci should have specific segregation patterns based on marker type.

# calculate genotype frequencies and p-value (after Bonferroni correction)
geno.freq.l <- geno.table(FamA_M.f3)

# create list of markers with p-value <0.05
seg.dist.l <- geno.freq.l[geno.freq.l$P.value < 0.05/totmar(FamA_M.f3),] 
drop.l.seg.dist <- rownames(geno.freq.l[geno.freq.l$P.value < 1e-10,])

# drop markers with significant (p<0.05) deviation from expected segregation patterns
FamA_M.f4 <- drop.markers(FamA_M.f3, drop.l.seg.dist)

# number of individuals/loci in data set
Sum_FamA_M.f4 <- summary(FamA_M.f4)
Nind_M.f4 <- as.numeric(Sum_FamA_M.f4$n.ind)
Nloci_M.f4 <- sum(as.numeric(Sum_FamA_M.f4$n.mar))

After filtering 1241 loci are retained.

Filter 5: Genotype frequencies by individual

Determine genotype frequencies by individual and remove individuals as necessary.

# if no individuals removed
FamA_M.f5 <- FamA_M.f4

# number of individuals/loci in data set
Sum_FamA_M.f5 <- summary(FamA_M.f5)
Nind_M.f5 <- as.numeric(Sum_FamA_M.f5$n.ind)
Nloci_M.f5 <- sum(as.numeric(Sum_FamA_M.f5$n.mar))

After filtering 152 individuals are retained.

Determine correct linkage phase

Estimate recombination fraction (rf) & LOD score for each marker pair. Potentially switched markers are pairs of marker with strong association (LOD score) but rf >> 1/4.

Form initial linkage groups and re-organize loci into initial LGs.

# form initial linkage groups
init.lg <- formLinkageGroups(FamA_M.LP, max.rf = 0.35, min.lod = 6)
table(init.lg[,2])

 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 
30 29 28 25 25 25 25 24 24 23 23 22 22 21 21 21 21 21 20 19 19 19 19 18 18 18 18 18 
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 
17 17 16 16 16 16 16 16 15 15 14 14 14 14 13 13 13 12 12 12  3  1  1 

# re-organize markers into inferred initial LGs
FamA_M.LP <- formLinkageGroups(FamA_M.LP, max.rf = 0.35, min.lod = 6, reorgMarkers=TRUE)

Plot rf vs LOD to identify chromosomes to switch alleles for.

Switch alleles, re-estimate rf, form new linkage groups and re-organize markers.

# use plot to identify chromosomes to switch alleles for
for (chr in c(6, 12:14, 16:19,
              26, 28:30, 35, 37:41, 43:48)){
  toswitch <- markernames(FamA_M.LP, chr=(chr))
  FamA_M.LP <- switchAlleles(FamA_M.LP, toswitch)
}

# re-estimate rf
FamA_M.LP <- est.rf(FamA_M.LP)

# form new linkage groups & re-organize markers
init.lg <- formLinkageGroups(FamA_M.LP, max.rf=0.35, min.lod=6)
table(init.lg[,2])
FamA_M.LP <- formLinkageGroups(FamA_M.LP, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)

Plot rf vs LOD after switching alleles:

# plot rf vs LOD to see re-organized markers
heatMap(FamA_M.LP, alternate.chrid = TRUE)

Switch alleles, re-estimate rf, form new linkage groups and re-organize markers.

# use plot to identify chromosomes to switch alleles for
for (chr in c(24, 26)){
  toswitch <- markernames(FamA_M.LP, chr=(chr))
  FamA_M.LP <- switchAlleles(FamA_M.LP, toswitch)
}

# re-estimate rf
FamA_M.LP <- est.rf(FamA_M.LP)

# form new linkage groups & re-organize markers
init.lg <- formLinkageGroups(FamA_M.LP, max.rf=0.35, min.lod=6)
table(init.lg[,2])
FamA_M.LP <- formLinkageGroups(FamA_M.LP, max.rf=0.35, min.lod=6, reorgMarkers=TRUE)

Plot rf vs LOD after switching alleles:

Drop marker in LG25 and 26:

pull.map(FamA_M.LP, chr="26")

dDocent_Contig_25489 
                   0 

Establish initial marker order

Establish initial marker order within LGs using greedy algorithm. Adds one marker at a time with previously added marker locations fixed. Position determined by minimizing number of obligate cross overs.

# Create final cross/map object
FamA_M.final <- FamA_M.LP

# establish initial marker orders
FamA_M.final <- orderMarkers(FamA_M.final, chr = 1:24)

Initial quality assessment

Examine pairwise linkage & Identify gaps between markers

Plot rf vs LOD after ordering markers:

Large gaps between markers indicate that adjacent markers are only weakly linked. Compare average and largest spacing per LG.

summaryMap(FamA_M.final)

        n.mar length ave.spacing max.spacing
1          53   68.9         1.3         6.6
2          50   67.5         1.4         7.3
3          49  154.2         3.2        78.6
4          46   57.6         1.3        14.2
5          46  104.8         2.3        20.3
6          45   86.9         2.0         9.5
7          44  119.8         2.8        18.5
8          43   69.4         1.7        12.6
9          42   56.8         1.4         8.5
10         41   78.3         2.0        12.6
11         39   82.3         2.2        10.2
12         37   48.7         1.4         7.3
13         36   87.6         2.5        17.6
14         33   65.7         2.1        17.6
15         33   66.7         2.1        11.8
16         33   50.4         1.6         7.0
17         32   59.4         1.9        15.0
18         32   51.5         1.7         8.7
19         31  100.9         3.4        36.5
20         31   63.5         2.1        11.0
21         29   55.9         2.0         8.7
22         29   58.6         2.1        18.5
23         28   94.5         3.5        43.4
24         28   74.3         2.8        17.6
overall   910 1824.3         2.1        78.6

Drop markers one by one to determine change in chromosome length & change in log likelihood. Identify markers that result in smaller estimated chromosome length and increased LOD. Dropping terminal markers will usually result in decreased length.

# drop each marker and compare effect
dropone <- droponemarker(FamA_M.final, maxit = 4000, tol = 1e-6, 
                         error.prob = 0.001, verbose = TRUE) %>%
  add_rownames("Marker") %>%
  gather("Effect", "Change", 4:5)

# plot results
plot_dropone <- function(dropone){
  p <- ggplot(dropone, aes(x = pos, y = Change)) +
    geom_line() +
    geom_point(shape = 21, size = 2,
               color = "black", fill = "white") +
    facet_grid(Effect ~ ., scales = "free") +
    labs(x = paste("LG", dropone$chr, "Position", sep =" "), y = "Change if marker dropped") +
    theme_facet
  print(p)
}

dropone %>%
  group_by(chr) %>%
  do(plot = plot_dropone(.)) %>%
  ungroup()

Compare effect in terms of change in length and likelihood and identify top 10 loci with strongest effect per LG.

# data frame with top10 markers with strongest effect
dropone_sum <- dropone %>%
  spread("Effect", "Change") %>%
  group_by(chr) %>%
  top_n(10, LOD) %>%
  ungroup() %>%
  arrange(chr, LOD)

Drop markers as necessary and re-estimate map.

Identify individuals with high number of XO

Observed number of crossovers per individual; eliminate individuals with excessively high number of crossovers >50).

Estimate genotyping error rate

Estimate genotyping error rate.

# Estimate genotyping error rate
loglik <- err <- c(0.0001, 0.0005, 0.001, 0.0025, 0.005, 0.0075, 0.01, 0.015, 0.02)
for(i in seq(along=err)) {
  cat(i, "of", length(err), "\n")
  tempmap <- est.map(FamA_M.final, error.prob=err[i])
  loglik[i] <- sum(sapply(tempmap, attr, "loglik"))
}

1 of 9 
2 of 9 
3 of 9 
4 of 9 
5 of 9 
6 of 9 
7 of 9 
8 of 9 
9 of 9 

lod <- (loglik - max(loglik))/log(10)
# plot log10 likelihood
par(mfrow=c(1,1), ps=12, cex=1)
plot(err, lod, xlab="Genotyping error rate", xlim=c(0,0.02), ylab=expression(paste(log[10], "likelihood")))

Genotyping error rate with highest likelihood is 0.0005.

Identify genotyping errors and code as missing

Tight double-crossovers (single marker out of phase with adjacent markers) indicate genotyping errors. Identify double cross-overs by calculating genotyping error LOD scores (Lincoln & Lander 1992). Calculate errorload (genotype being in error vs not being in error). Plot grid of LOD scores indicating genotypes likely to be in error: darker pixels likely to be error: white (LOD <2), grey (2 < LOD < 3), pink (3 < LOD < 4.5), purple (LOD >4.5)

Eliminate genotypes (coded as missing) for LOD > 6 and re-estimate map.

Re-order markers by LG and compare alternate orders

Compare number obligate XO, LOD & Chromosome length to determine best order

Re-order markers after loci and individuals have been removed and use sliding window to compare alternate orders in terms of number of obligate cross overs and likelihood.

# compare number of obligate crossovers for alternate orders
FamA_M.XO <- lapply(c(1:24), function(chr) {
  rip <- ripple(FamA_M.final, chr = chr, window = 7, method="countxo", error.prob = 0.001, maxit = 2000, tol = 1e-6)
  rip <- rip[1:11,]
})

# access each element in list using [[]]
# Ch1_XO <- as.data.frame(FamA_M.XO[[1]]) %>%
#  rename(Chr1 = obligXO) %>%
#  select(Chr1)

FamA_M.lik <- lapply(c(1:24), function(chr) {
  rip.lik <- ripple(FamA_M.final, chr = chr, window = 5, method="likelihood", error.prob=0.002, maxit=2000, tol=1e-4)
  rip.lik <- rip.lik[1:11,]
})

Switch marker order as necessary

Identify order with lowest number of XO and highest LOD while minimizing chromosome length.

change.marker.order <- function(LG, neworder) {
  Orderlikelih <- as.data.frame(FamA_M.lik[[LG]]) %>%
    select(-LOD) %>%
    select(-chrlen)
  switch_marker <- as.numeric(Orderlikelih[3,])
  FamA_M.final <- switch.order(FamA_M.final, chr = LG, switch_marker, error.prob=0.001)
}

## LG1
# Identify order with lowest no. of XO & highest LOD while minimizing chromosome length
FamA_M.XO[[1]]
FamA_M.lik[[1]]
# change.marker.order(1, new)

## LG2 
FamA_M.XO[[2]]
FamA_M.lik[[2]]
# change.marker.order(2, 1)

## LG3
FamA_M.XO[[3]]
FamA_M.lik[[3]]
# change.marker.order(3, 1)

## LG4
FamA_M.XO[[4]]
FamA_M.lik[[4]]
# change.marker.order(4, neworder)

## LG5
FamA_M.XO[[5]]
FamA_M.lik[[5]]
# change.marker.order(5, 2)

## LG6
FamA_M.XO[[6]]
FamA_M.lik[[6]]
# change.marker.order(6, neworder)

## LG7 
FamA_M.XO[[7]]
FamA_M.lik[[7]]
# change.marker.order(7, neworder)

## LG8
FamA_M.XO[[8]]
FamA_M.lik[[8]]
# change.marker.order(8, neworder)

## LG9
FamA_M.XO[[9]]
FamA_M.lik[[9]]
# change.marker.order(9, neworder)

## LG10
FamA_M.XO[[10]]
FamA_M.lik[[10]]
# change.marker.order(10, 1)

## LG11
FamA_M.XO[[11]]
FamA_M.lik[[11]]
change.marker.order(11, 1)

## LG12
FamA_M.XO[[12]]
FamA_M.lik[[12]]
# change.marker.order(12, 1)

## LG13
FamA_M.XO[[13]]
FamA_M.lik[[13]]
# change.marker.order(13, new)

## LG14
FamA_M.XO[[14]]
FamA_M.lik[[14]]
# change.marker.order(14, new)

## LG15
FamA_M.XO[[15]]
FamA_M.lik[[15]]
# change.marker.order(15, 1)

## LG16
FamA_M.XO[[16]]
FamA_M.lik[[16]]
# change.marker.order(16, neworder)

## LG17
FamA_M.XO[[17]]
FamA_M.lik[[17]]
# change.marker.order(17, neworder)

## LG18
FamA_M.XO[[18]]
FamA_M.lik[[18]]
# change.marker.order(18, 1)

## LG19
FamA_M.XO[[19]]
FamA_M.lik[[19]]
# change.marker.order(19, new)

## LG20
FamA_M.XO[[20]]
FamA_M.lik[[20]]
# change.marker.order(20, new)

## LG21
FamA_M.XO[[21]]
FamA_M.lik[[21]]
# change.marker.order(21, 4)

## LG22
FamA_M.XO[[22]]
FamA_M.lik[[22]]
# change.marker.order(22, new)

## LG23
FamA_M.XO[[23]]
FamA_M.lik[[23]]
# change.marker.order(23, new)

## LG24
FamA_M.XO[[24]]
FamA_M.lik[[24]]
# change.marker.order(24, new)

# re-estimate map after re-ordering all LGs (chromosomes)
FamA_M.reorder <- est.map(FamA_M.final, error.prob = 0.001)

# re-place map portion of cross
FamA_M.reorder <- replace.map(FamA_M.final, FamA_M.reorder)

Post-odering quality assessment II

Examine pairwise linkage & Identify gaps between markers

Plot rf vs LOD

Identify large gaps between markers which is indicative of adjacent markers being only weakly linked.

summaryMap(FamA_M.final)

        n.mar length ave.spacing max.spacing
1          52   63.0         1.2         6.6
2          49   61.7         1.3         7.3
3          47   62.0         1.3        11.8
4          46   57.6         1.3        14.2
5          45   92.1         2.1        20.3
6          44   84.4         2.0         9.5
7          43  119.2         2.8        29.4
8          43   64.2         1.5        12.6
9          42   56.7         1.4         8.4
10         41   77.6         1.9        12.6
11         38   73.5         2.0        10.2
12         37   48.6         1.4         7.3
13         35   93.8         2.8        36.6
14         33   63.6         2.0        17.6
15         32   64.7         2.1        15.9
16         33   49.0         1.5         7.0
17         31   49.7         1.7        14.8
18         32   50.4         1.6         8.8
19         30   64.3         2.2        13.4
20         31   62.3         2.1        10.8
21         29   55.2         2.0         8.7
22         28   56.6         2.1        18.5
23         27   49.9         1.9         8.4
24         28   72.8         2.7        17.6
overall   896 1593.0         1.8        36.6

Plot maps of all linkage groups:

Drop markers one by one determine change in chromosome length & change in log likelihood. Identify markers that result in smaller estimated chromosome length and increased LOD.

# drop each marker and compare effect
dropone <- droponemarker(FamA_M.final, maxit = 4000, tol = 1e-6, 
                         error.prob = 0.001, verbose = TRUE) %>%
  add_rownames("Marker") %>%
  gather("Effect", "Change", 4:5)

# plot results
plot_dropone <- function(dropone){
  p <- ggplot(dropone, aes(x = pos, y = Change)) +
    geom_line() +
    geom_point(shape = 21, size = 2,
               color = "black", fill = "white") +
    facet_grid(Effect ~ ., scales = "free") +
    labs(x = paste("LG", dropone$chr, "Position", sep =" "), y = "Change if marker dropped") +
    theme_facet
  print(p)
}

dropone %>%
  group_by(chr) %>%
  do(plot = plot_dropone(.)) %>%
  ungroup()

Compare effect in terms of change in length and likelihood and identify top 10 loci with strongest effect per LG.

# data frame with top10 markers with strongest effect
dropone_sum <- dropone %>%
  spread("Effect", "Change") %>%
  group_by(chr) %>%
  top_n(10, LOD) %>%
  ungroup() %>%
  arrange(chr, LOD)

Identify individuals with high number of XO

Finalize map

# print summary map and map data
SummaryFamA_M <- summaryMap(FamA_M.final)

write.table(SummaryFamA_M, file = "results/SummaryFamA_M.map",
          quote = FALSE, sep = "\t", row.names = FALSE)

FamA_M.map <- pull.map(FamA_M.final,as.table = TRUE)

write.table(FamA_M.map, file = "data/LINKMAP/FamA_M.map",
            quote = FALSE, sep = "\t", row.names = TRUE)

Write files of finalized map. Map contains 1311 loci grouped into 24 LG.

Compare male/female maps

Re-insert matching loci (same segregation pattern)

Finalize the female map by inserting the loci that had duplicate segregation patterns that were previously removed.

# import female map
FamA_F.map <- read.table("data/LINKMAP/FamA_F.map", stringsAsFactors = FALSE, header = FALSE,
                         col.names = c("LOCUS", "LG_F", "POS_F"), skip = 1)
# add duplicate loci back in
match.loci_F <- read.table("data/LINKMAP/MatchLociFamAF.txt", header = TRUE,
                         stringsAsFactors = FALSE) %>%
  left_join(FamA_F.map) %>%
  select(MATCHLOC, LG_F, POS_F) %>%
  rename(LOCUS = MATCHLOC)

Joining, by = "LOCUS"

# update female map
FamA_F.map <- bind_rows(FamA_F.map, match.loci_F) %>%
  group_by(LG_F) %>%
  arrange(POS_F)
write.table(FamA_F.map, file = "results/FamA_F.map",
          quote = FALSE, sep = "\t", row.names = FALSE)

Finalize the male map by inserting the loci that had duplicate segregation patterns that were previously removed.

# import male map
FamA_M.map <- read.table("data/LINKMAP/FamA_M.map", stringsAsFactors = FALSE, header = FALSE,
                         col.names = c("LOCUS", "LG_M", "POS_M"), skip = 1)
# add duplicate loci back in
match.loci_M <- read.table("data/LINKMAP/MatchLociFamAM.txt", header = TRUE,
                         stringsAsFactors = FALSE) %>%
  left_join(FamA_M.map) %>%
  select(MATCHLOC, LG_M, POS_M) %>%
  rename(LOCUS = MATCHLOC)

Joining, by = "LOCUS"

# update female map
FamA_M.map <- bind_rows(FamA_M.map, match.loci_M) %>%
  group_by(LG_M) %>%
  arrange(POS_M)
write.table(FamA_M.map, file = "results/FamA_M.map",
          quote = FALSE, sep = "\t", row.names = FALSE)

Compare & join male/female maps

M.unique

[1] 477

For mapping family A, the female and male maps contain 1321 and 1311 loci, respectively. Of those the 834 are shared, while the female and male maps contain 487 and 477 unique loci, respectively.

FamA family map (onemap)

Map all loci segregating in male and/or female using full-cross.

Filter onemap input file by markers in joined map data set

Import tab delimited onemap input file converted from filtered *.tsv-file using map_convert.pl.

# change header names to Marker | Markertype | Genotype
FamA.onemap <- read.delim("data/LINKMAP/FamA.onemap", header = FALSE,
                          col.names = c("LOCUS", "LOCUSTYPE", "GENOTYPE"),
                          skip = 1, stringsAsFactors = FALSE)
# number of individuals in dataset
N <- read.delim("data/LINKMAP/FamA.onemap", header = FALSE,
                          col.names = c("N_IND", "N_LOCI"),
                          nrows = 1, stringsAsFactors = FALSE)
# keep only markers in combined data set (filter by joined data set)
Onemap.filtered <- semi_join(FamA.onemap, mapped_markers)

Joining, by = "LOCUS"

Unfiltered onemap-data set contains 162 individuals genotyped at 1939 loci. After removing loci filtered during creation of male and female map, the data set consists of nrow(Onemap.filtered) loci.

Markers with matching genotypes (same segregation pattern) are removed and stored in table with corresponding loci retained in data set to insert after mapping loci.

# remove duplicate genotypes
Onemap.filt2 <- Onemap.filtered %>%
  distinct(GENOTYPE, .keep_all = TRUE)
# create dataframe of removed loci and corresponding loci
match_loc <- anti_join(Onemap.filtered, Onemap.filt2) %>%
  rename(MATCHLOC = LOCUS) %>%
  select(MATCHLOC, GENOTYPE) %>% 
  left_join(Onemap.filt2, by = "GENOTYPE") %>%
  select(LOCUS, MATCHLOC)

Joining, by = c("LOCUS", "LOCUSTYPE", "GENOTYPE")

# write filtered input file for Onemap
# need to add extra line back in
con <- file("data/LINKMAP/FamAfilter.onemap", open = "wt")
writeLines(paste(N$N_IND, (nrow(Onemap.filt2))), con)
write.table(Onemap.filt2, con, 
            append = FALSE, quote = FALSE, sep = " ", 
            na = "NA", dec = ".", 
            col.names = FALSE, row.names = FALSE)
close(con)

The filtered data set contains 1565 loci.

Use onemap to create linkage map based on full-siblings

Before importing use information from *xo file to remove genotypes flagged as potential genotyping error in filtered onemap input.

FamA.Map

This is an object of class 'outcross'
    No. individuals:    162 
    No. markers:        1565 
    Segregation types:
       A.1: 389
       A.2: 483
       D1.10:   346
       D2.15:   347

Estimate two-point recombination fractions between all pairs of markers (Wu et al. 2002).

FamA.Map.rf.L6.rf35 <- rf.2pts(FamA.Map, LOD = 6, max.rf = 0.35, verbose=TRUE)

  |                                                                                         
  |                                                                                   |   0%
  |                                                                                         
  |                                                                                   |   1%
  |                                                                                         
  |=                                                                                  |   1%
  |                                                                                         
  |=                                                                                  |   2%
  |                                                                                         
  |==                                                                                 |   2%
  |                                                                                         
  |==                                                                                 |   3%
  |                                                                                         
  |===                                                                                |   3%
  |                                                                                         
  |===                                                                                |   4%
  |                                                                                         
  |====                                                                               |   4%
  |                                                                                         
  |====                                                                               |   5%
  |                                                                                         
  |=====                                                                              |   5%
  |                                                                                         
  |=====                                                                              |   6%
  |                                                                                         
  |=====                                                                              |   7%
  |                                                                                         
  |======                                                                             |   7%
  |                                                                                         
  |======                                                                             |   8%
  |                                                                                         
  |=======                                                                            |   8%
  |                                                                                         
  |=======                                                                            |   9%
  |                                                                                         
  |========                                                                           |   9%
  |                                                                                         
  |========                                                                           |  10%
  |                                                                                         
  |=========                                                                          |  10%
  |                                                                                         
  |=========                                                                          |  11%
  |                                                                                         
  |==========                                                                         |  11%
  |                                                                                         
  |==========                                                                         |  12%
  |                                                                                         
  |==========                                                                         |  13%
  |                                                                                         
  |===========                                                                        |  13%
  |                                                                                         
  |===========                                                                        |  14%
  |                                                                                         
  |============                                                                       |  14%
  |                                                                                         
  |============                                                                       |  15%
  |                                                                                         
  |=============                                                                      |  15%
  |                                                                                         
  |=============                                                                      |  16%
  |                                                                                         
  |==============                                                                     |  16%
  |                                                                                         
  |==============                                                                     |  17%
  |                                                                                         
  |===============                                                                    |  17%
  |                                                                                         
  |===============                                                                    |  18%
  |                                                                                         
  |===============                                                                    |  19%
  |                                                                                         
  |================                                                                   |  19%
  |                                                                                         
  |================                                                                   |  20%
  |                                                                                         
  |=================                                                                  |  20%
  |                                                                                         
  |=================                                                                  |  21%
  |                                                                                         
  |==================                                                                 |  21%
  |                                                                                         
  |==================                                                                 |  22%
  |                                                                                         
  |===================                                                                |  22%
  |                                                                                         
  |===================                                                                |  23%
  |                                                                                         
  |====================                                                               |  23%
  |                                                                                         
  |====================                                                               |  24%
  |                                                                                         
  |====================                                                               |  25%
  |                                                                                         
  |=====================                                                              |  25%
  |                                                                                         
  |=====================                                                              |  26%
  |                                                                                         
  |======================                                                             |  26%
  |                                                                                         
  |======================                                                             |  27%
  |                                                                                         
  |=======================                                                            |  27%
  |                                                                                         
  |=======================                                                            |  28%
  |                                                                                         
  |========================                                                           |  28%
  |                                                                                         
  |========================                                                           |  29%
  |                                                                                         
  |========================                                                           |  30%
  |                                                                                         
  |=========================                                                          |  30%
  |                                                                                         
  |=========================                                                          |  31%
  |                                                                                         
  |==========================                                                         |  31%
  |                                                                                         
  |==========================                                                         |  32%
  |                                                                                         
  |===========================                                                        |  32%
  |                                                                                         
  |===========================                                                        |  33%
  |                                                                                         
  |============================                                                       |  33%
  |                                                                                         
  |============================                                                       |  34%
  |                                                                                         
  |=============================                                                      |  34%
  |                                                                                         
  |=============================                                                      |  35%
  |                                                                                         
  |=============================                                                      |  36%
  |                                                                                         
  |==============================                                                     |  36%
  |                                                                                         
  |==============================                                                     |  37%
  |                                                                                         
  |===============================                                                    |  37%
  |                                                                                         
  |===============================                                                    |  38%
  |                                                                                         
  |================================                                                   |  38%
  |                                                                                         
  |================================                                                   |  39%
  |                                                                                         
  |=================================                                                  |  39%
  |                                                                                         
  |=================================                                                  |  40%
  |                                                                                         
  |==================================                                                 |  40%
  |                                                                                         
  |==================================                                                 |  41%
  |                                                                                         
  |==================================                                                 |  42%
  |                                                                                         
  |===================================                                                |  42%
  |                                                                                         
  |===================================                                                |  43%
  |                                                                                         
  |====================================                                               |  43%
  |                                                                                         
  |====================================                                               |  44%
  |                                                                                         
  |=====================================                                              |  44%
  |                                                                                         
  |=====================================                                              |  45%
  |                                                                                         
  |======================================                                             |  45%
  |                                                                                         
  |======================================                                             |  46%
  |                                                                                         
  |=======================================                                            |  46%
  |                                                                                         
  |=======================================                                            |  47%
  |                                                                                         
  |=======================================                                            |  48%
  |                                                                                         
  |========================================                                           |  48%
  |                                                                                         
  |========================================                                           |  49%
  |                                                                                         
  |=========================================                                          |  49%
  |                                                                                         
  |=========================================                                          |  50%
  |                                                                                         
  |==========================================                                         |  50%
  |                                                                                         
  |==========================================                                         |  51%
  |                                                                                         
  |===========================================                                        |  51%
  |                                                                                         
  |===========================================                                        |  52%
  |                                                                                         
  |============================================                                       |  52%
  |                                                                                         
  |============================================                                       |  53%
  |                                                                                         
  |============================================                                       |  54%
  |                                                                                         
  |=============================================                                      |  54%
  |                                                                                         
  |=============================================                                      |  55%
  |                                                                                         
  |==============================================                                     |  55%
  |                                                                                         
  |==============================================                                     |  56%
  |                                                                                         
  |===============================================                                    |  56%
  |                                                                                         
  |===============================================                                    |  57%
  |                                                                                         
  |================================================                                   |  57%
  |                                                                                         
  |================================================                                   |  58%
  |                                                                                         
  |=================================================                                  |  58%
  |                                                                                         
  |=================================================                                  |  59%
  |                                                                                         
  |=================================================                                  |  60%
  |                                                                                         
  |==================================================                                 |  60%
  |                                                                                         
  |==================================================                                 |  61%
  |                                                                                         
  |===================================================                                |  61%
  |                                                                                         
  |===================================================                                |  62%
  |                                                                                         
  |====================================================                               |  62%
  |                                                                                         
  |====================================================                               |  63%
  |                                                                                         
  |=====================================================                              |  63%
  |                                                                                         
  |=====================================================                              |  64%
  |                                                                                         
  |======================================================                             |  64%
  |                                                                                         
  |======================================================                             |  65%
  |                                                                                         
  |======================================================                             |  66%
  |                                                                                         
  |=======================================================                            |  66%
  |                                                                                         
  |=======================================================                            |  67%
  |                                                                                         
  |========================================================                           |  67%
  |                                                                                         
  |========================================================                           |  68%
  |                                                                                         
  |=========================================================                          |  68%
  |                                                                                         
  |=========================================================                          |  69%
  |                                                                                         
  |==========================================================                         |  69%
  |                                                                                         
  |==========================================================                         |  70%
  |                                                                                         
  |===========================================================                        |  70%
  |                                                                                         
  |===========================================================                        |  71%
  |                                                                                         
  |===========================================================                        |  72%
  |                                                                                         
  |============================================================                       |  72%
  |                                                                                         
  |============================================================                       |  73%
  |                                                                                         
  |=============================================================                      |  73%
  |                                                                                         
  |=============================================================                      |  74%
  |                                                                                         
  |==============================================================                     |  74%
  |                                                                                         
  |==============================================================                     |  75%
  |                                                                                         
  |===============================================================                    |  75%
  |                                                                                         
  |===============================================================                    |  76%
  |                                                                                         
  |===============================================================                    |  77%
  |                                                                                         
  |================================================================                   |  77%
  |                                                                                         
  |================================================================                   |  78%
  |                                                                                         
  |=================================================================                  |  78%
  |                                                                                         
  |=================================================================                  |  79%
  |                                                                                         
  |==================================================================                 |  79%
  |                                                                                         
  |==================================================================                 |  80%
  |                                                                                         
  |===================================================================                |  80%
  |                                                                                         
  |===================================================================                |  81%
  |                                                                                         
  |====================================================================               |  81%
  |                                                                                         
  |====================================================================               |  82%
  |                                                                                         
  |====================================================================               |  83%
  |                                                                                         
  |=====================================================================              |  83%
  |                                                                                         
  |=====================================================================              |  84%
  |                                                                                         
  |======================================================================             |  84%
  |                                                                                         
  |======================================================================             |  85%
  |                                                                                         
  |=======================================================================            |  85%
  |                                                                                         
  |=======================================================================            |  86%
  |                                                                                         
  |========================================================================           |  86%
  |                                                                                         
  |========================================================================           |  87%
  |                                                                                         
  |=========================================================================          |  87%
  |                                                                                         
  |=========================================================================          |  88%
  |                                                                                         
  |=========================================================================          |  89%
  |                                                                                         
  |==========================================================================         |  89%
  |                                                                                         
  |==========================================================================         |  90%
  |                                                                                         
  |===========================================================================        |  90%
  |                                                                                         
  |===========================================================================        |  91%
  |                                                                                         
  |============================================================================       |  91%
  |                                                                                         
  |============================================================================       |  92%
  |                                                                                         
  |=============================================================================      |  92%
  |                                                                                         
  |=============================================================================      |  93%
  |                                                                                         
  |==============================================================================     |  93%
  |                                                                                         
  |==============================================================================     |  94%
  |                                                                                         
  |==============================================================================     |  95%
  |                                                                                         
  |===============================================================================    |  95%
  |                                                                                         
  |===============================================================================    |  96%
  |                                                                                         
  |================================================================================   |  96%
  |                                                                                         
  |================================================================================   |  97%
  |                                                                                         
  |=================================================================================  |  97%
  |                                                                                         
  |=================================================================================  |  98%
  |                                                                                         
  |================================================================================== |  98%
  |                                                                                         
  |================================================================================== |  99%
  |                                                                                         
  |===================================================================================|  99%
  |                                                                                         
  |===================================================================================| 100%

Assign markers to linkage groups. Create sequence of all markers to assign to LGs and assign LG for LOD > 6 and rf < 0.35.

Markers are grouped into 24 LGs.

Map markers within linkage group using haldane mapping function. Onemap first creates a framework consisting of the most informative markers (segregating in both parents) and then maps the remaining markers.

# define map function
set.map.fun(typ = "haldane")

# create list of separate LGs 
# (convert from class group to class sequence)
LG_list <- lapply(c(1:24), function(LG) {
  make.seq(LGs, LG)
})

# uses compare() to create a framework of most informative markers
# map remaining markers using try.seq()
# this step will take a few hours
LG_ordered <- lapply((LG_list), function(LG){
  order.seq(LG, n.init=6, THRES=3, touchdown=TRUE, tol=10e-2)
})

# create marker order with all markers
LG_final <- lapply((LG_ordered), function(LG){
  make.seq(LG, "force")
})

Compare alternative marker sequences

Compare alternative sequences using ripple for 4 marker window. Make changes as needed (more positive LOD indicates better sequence).

# compare and display alternative orders for LG1
# prints result of each window to screen
# more positive LOD indicates better sequence
# set window size using ws, LOD sets threshold for alt orders printed

# Start writing to an output file
sink("data/LINKMAP/ripple.onemap")

# run ripple for each LG
lapply((LG_final), function(LG){
  ripple.seq(LG, ws = 4, LOD = 3, tol = 10e-2)
})

# Stop writing to the file
sink()

# change marker order if necessary for each LG# using 
# LG#.new.seq <- make.seq(LGx.final, c())
# LG1.final <- map(LG#.new.seq)

Plot recombination fraction matrix to inspect linkage between mapped markers.

Draw linkage map and save results

Add duplicate markers (with same genotype/segregation pattern) and save finalized map to file